! ==============================================================================
! This file is part of CUPP (Community Unified Physics Package) since 2023.
! ==============================================================================

module cupp_math_mod

  use cupp_const_mod

  implicit none

contains

  pure elemental real(r8) function f_gamma(x) result(res)

    real(r8), intent(in) :: x

    real(r8), parameter :: b(8) = [                                  &
      -0.577191652_r8,0.988205891_r8,-0.897056937_r8,0.918206857_r8, &
      -0.756704078_r8,0.482199394_r8,-0.193527818_r8,0.035868343_r8  &
    ]
    real(r8) f, g, y
    integer i

    f = 1
    y = x
    do i = 1, 200
      if (y <= 2) exit
      y = y - 1
      f = f * y
    end do
    g = 1
    y = y - 1
    do i = 1, size(b)
      g = g + b(i) * y**i
    end do
    res = f * g

  end function f_gamma

end module cupp_math_mod